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Abstract. We describe the numerical code N-MODY, a parallel particle-mesh code for 
collisionless N-body simulations in modified Newtonian dynamics (MOND). N-MODY is 
based on a numerical potential solver in spherical coordinates that solves the non-linear 
MOND field equation, and is ideally suited to simulate isolated stellar systems. N-MODY 
can be used also to compute the MOND potential of arbitrary static density distributions. 
A few applications of N-MODY indicate that some astrophysically relevant dynamical pro- 
cesses are profoundly different in MOND and in Newtonian gravity with dark matter. 
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1. Introduction 

Modified Ne wtonia n dynamics (MOND) is an alternative gravity theory, originally proposed 
by iMilgrom] d 1983b to explain the observed kinematics of dis k galaxies without dark matte r. In 
Bekenstein & Milgrom's Lagrangian formulation of MOND (Bekenstein & Milgrom 1984) the 
Poisson equation is replaced by the non-linear field equation 



where «o - 1.2 x 10~ I0 ms~ 2 is a characteristic acceleration, ||...|| is the standard Euclidean norm 
in R 3 , and (f> is the MOND gravitational potential produced by the density distribution p. The 
MOND interpolating function //(y) is a monotonic function such that 
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For finite mass systems the boundary condition of equation (HJ is — » for ||x|| — > oo. 

As MOND i s successful in reproducing several observed properties of galaxies (e.g. 
Bek ensteinl I200 6). there is growing interest in studying dynamical processes in MOND with 
the aid of N-body simulations. While in Newtonian gravity N-body codes can take advantage o f 
the multipole expansion of the Green function (treecodes: iBarnes & Huj[l986t lDehnenll2002l) . 
the non-linearity of the MOND fiel d equation forces one to solv e for the potential on a grid, as 
done in particle-mesh schemes (see lHockney & Eastwood 1981). 

We developed N-MODY, a parallel three-dimensional particle-mesh code for collisionless N- 
body simulations in MOND. The N-body code and the poten t ial sol ver o n which the c o de is ba sed 
have been already tested and applied in ICiotti et alj d2006l [20071) and INipoti et alj ([2007a b c, 
2008). The potential solver of N-MODY solves the MOND field equation O using a relaxation 
method in spherical coordinates based on the spherical harmoni cs expansion. Thus the code is 
ideally suited for simulations of isolated stellar systems (but see INipoti et al]|2007cl for an ap- 
plication to galaxy merging). N-MODY can also be used to compute the MOND potential of 
arbitrary static density distributions. N-MODY is one of the very few MOND N-body code de- 
veloped so far: as far as we know the only othe r three-dimensional MOND N-body code is Brada 
and Milgrom's code (Brada & Milgrom 1999), which is based on a multi-grid potential so lver in 
Cartesian coordinates and has been recently implemented also bv lTiret & Combesl d2007l) . 

2. Overview of the code 

N-MODY implements a particle-mesh scheme in spherical coordinates, following these main 
computational steps: 

1) for a given distribution of N particles, a grid-based density field is reconstructed by mass 
deposition with linear or quadratic shape functions. Particles are represented by a six- 
component array (x p , \ p ) of positions and velocities in Cartesian components (p — 1 , ...,N). 

2) For a given density distribution, the MOND acceleration and potential fields are computed 
on the grid points. As an alternative, the code also provides a fast Newtonian solver. 

3) To move particles, the spherical components of the acceleration are first interpolated at each 
particle position (using the same linear or quadratic shape functions used in the mass deposi- 
tion) and then transformed into Cartesian components. 

4) Finally, particle positions and velocities are advanced in time using a leapfrog scheme of 
either second or fourth order. 

In N-MODY steps 1) and 3) are parallelized using MPI routines: particles are uniformly 
distributed among the processors (PEs) following the sequential ordering provided by their initial 
memory addresses and then never exchanged among different PEs. This simple strategy assures 
a full efficiency in parallel execution, but entails more memory resources than in a standard 
domain decomposition. In fact, grid data computed on the overall computational domain must be 
at disposal of each PE at each timestep for particles interpolation and move. Step 2), which is only 
partially parallelized, contains the MOND potential solver, which represents a new contribution 
to the fast numerical solution of a non-linear elliptic equation. We now discuss in more detail the 
different parts of the code: the potential solver is described in Section[3] while the particle-mesh 
scheme and the time integration are described in Section[4] 

3. The MOND potential solver 

N-MODY solves the non-relativistic MOND field equation dT). By default in N-MODY we adopt 
as interpolating function p.{y) = yj sjl +y 2 , but it is clear that a trivial modification of the code 
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allows to implement any other continuous function /i with the required asymptotic behaviour^ 
N-MODY can also be used to simulate a Newtonian system, in which case the Poisson equation 
is solved instead of equation O, or a system in the so-called 'deep MOND regime', i.e. obeying 
the equation V • (||V0||V0) = 4nGa p. 

We consider only systems of finite mass M — J p(x) d 3 x, for which the boundary conditions 
of equation (UJ are given by ||V0|| — > 0(l/r) for r = \\x\\ — > oo. It must be stressed that the 
asymptotic behaviour of the MOND potential (<p ~ In r — > oo for r — » oo) is profoundly different 
from that of the Newtonian potential ((/>—> for r — » oo). For this reason, to solve numerically 
the MOND field equation (i) we discretize a sufficiently large computational domain, in a way 
the asymptotic boundary condition can be represented with reasonable accuracy and (ii) we use 
a relaxation method to solve the non-linear elliptic equation ((TJ with guess solution having the 
correct asymptotic behaviour. 



3.1. The computational grid 

To accomplish task (i) above, N-MODY uses a spherical grid (r, ■&, <p) with radial coordinate 
represented by the invertible mapping 



r®=Ltan a e, r'(# 



ah tan° 



cos 2 £ 



(3) 



where < £ < nJ2, the mapping index a - 1 or a — 2 and the scale length L are user provided 
parameters (see iLondrillo & MessinalH990l) . In this representation, the unbounded radial range 
(0, oo) is mapped onto the finite open interval [0, tt/2). The radial derivative is then expressed as 



d_ _ 1 d 

dr ~ r'(£) d% 

and the £ coordinate is discretized into the uniform grid 



(4) 



£ = (i + 1 /2)A£ A£ = — , i = 0, 1, ... N r , 



(5) 



so the corresponding discretized radial variable r, = r(^,) avoids the singular points r — and 
r = oo. The other coordinates <p) are discretized in the uniform grids 



#j = (j + 1 /2)A#, A# = n/N», j = Q,l,..,N*-l, 
to avoid the singular points = and ft = n, and 



(6) 



fe = 0, 1,..,A^-1 



(7) 



The corresponding , ^) numerical derivatives are computed by second-order or fourth-order 
finite difference central schemes. 



1 In versions of MOND based on Bekenstein's covariant theory TeVeS Bekensteii] 2004 ). there is a 
scalar field </> s obeying equation (TJ, but with an unbounded interpolating function n s instead of the bounded 
function ji. We verified that our code can be adapted to solve for S (see Famaev et alj2007 ). 
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3.2. A Newton-like relaxation procedure 

In N-MODY we use a Newton-like relaxation procedure to solve the non-linear MOND equation, 
which can be written as 



M[cf>(x)] = V • 



H - V0(x) 
fl() 



- AnGpix) = 0, g = 0(r _1 ) for r -> oo, 



(8) 



where p(x) is assigned, g(x) = ||g(x)||, with g(x) = -V0(x), and x = (r,d-,ip). Starting with an 
approximate guess solution (O) having the required asymptotic behaviour, a sequence of linear 
problems 



R 



-M[0 (n) ], n = 0,l, 



(9) 



is solved for the increments §<p w = <n+1) - <n) , each <$> {n) provisional solution having the same 
asymptotic behaviour as the guess (fP\ Here the choice of the relaxation linear operator R ( "^ is 
based on the requirement that it must assure convergence in the maximum norm ||...||oo, i.e. 



6TL = prAHL * 



a("-D|| 



n = 1,2,... 



(10) 



and it must be easy to invert. In a classical Newton method, one would use ^ (n) = 6M^"\ where 
the linear operator (5M (n) is such that 



5M [n) [54> {n) \ = M [4> (n+l) \ - M [</> (n) ] + O [(d<f> (n> ) 2 ] . 
For the specific case of M defined by equation ||5}, 
6M (n) = fi (n) V 2 + 6M { "\ 
where 

5M ( " ) = (V/i W ) • V + V ■ [n' (n) g (n \g {n) ■ v)] 



(11) 



(12) 



(13) 



and a . Boundedness of the inverse of the operator 6M (n) assures quadratic 

conve rgence of the scheme for </> (0) sufficiently close to the sought solution (e.g. lStoer & Bulirschl 
1 1 980t >. Unfortunately, [6M in ^Y l is difficult to compute, so we discretize the simpler linear oper- 
ator 



R (n) = ojfi {n) W 2 , 



(14) 



where lo > 1 is an empirical relaxation parameter. As an approximation of the Newton relax- 
ation operator 6M, this choice assures a lower convergence rate, but has clear computational 
advantages. Using equation (0 and the identity 



M(cf> (n) ) = 5M ( "- l) [S^ n - l) \ + M [(f> { "- l) ] + O [W ( " _1) ) 2 ] , 
it follows that the condition for convergence (TTOl i requires 



(") 



,(«) 



< 1. 



(15) 



(16) 



Thus N-MODY solves the sequence of Poisson equations: 



n = 0A, 



(17) 
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with source term given by 
1 



LO/J 



(») 



M 



In spherical coordinates, the Laplacian operator has the form 
1 



V 2 = 



where 



La = 



dr \ dr 



+ L§ + L v 



1 d 



sini? d& 



d 

sin § — 



i d 2 



sin 



2 §d>p 2 ' 



(18) 



(19) 



(20) 



After expanding the unknown function 6(p ( "\r, §, (p) and the source term «S (n) (r, tp) in spherical 
harmonics (or Fourier-Legendre) components 

5<p("\r, &,<p) = Y J W. <P)< S w (r, #.*») = 2 ^(r)Ff (-&, cp), (21) 

l.m l.m 

equation (fTTT i takes the simple one-dimensional form 
1 



4 p4i-^+d 



dr \ dr 



l.m 



(r) = rS { "'(r), 



(22) 



where we multiplied both sides by r to avoid the singularity in the source term for the astrophys- 
ically relevant case of p ~ r~ l central density profiles. Equation d22b . involving derivatives only 
in the radial coordinates, is discretized by central finite differences. 



3.3. Implementation of the relaxation procedure. 

The relaxation scheme in N-MODY is implemented in the following steps: 

1 . At the initial time t — in N-body simulations (and in the case of static models) the guess 
solution is chosen to be the spherically symmetric MOND solution of the angle-averaged 
density distribution 

Po oM = ~r I I P( r > <P) s i n d& dip, (23) 
4?rJo Jo 

which is easily derived from the corresponding spherical Newtonian solution 
dBekenstein & Milgroml 1 1 98-41) . In particular this guess solution satisfies the boundary 
condition, providing the values of the potential <f> <0> and of the radial acceleration at the 
boundary grid point r# r . At t > in N-body simulations the guess solution is provided by 
the numerical solution found in the previous timestep. 

2. For assigned acceleration g ( "'(r, §, tp) at the iteration level n, the source term rS ( "\r, §, <p) is 
evaluated, using central finite differences to approximate space derivatives. 

3. The source term is then transformed into Fourier-Legendre components by 

^(n) = ^ J e'^dip J SF>(r,, <p)P Um (fi) sin #d& = £ £ S { "\n, 0j, <pde- lm *P&(fij), (24) 
where PjUfij) is the inverse of the matrix Pi Jn (§j) of the associated Legendre polynomials, 
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4. The operator in equation d22l i is discretized in the radial coordinate by using finite differences 
for first and second derivatives. It results a tri-diagonal or penta-diagonal matrix of order N r + 
1 that can be easily inverted by using a standard Lower-Upper triangular (LU) decomposition 
to solve for the <WI variables with boundary conditions W„(Ov r ) = 0. 

5. The potential increments 8$ (r;) are back transformed into the (&, <p) coordinate space and 

the corresponding acceleration increments (5g (n) are evaluated using finite differences, to move 
to the next level solution g ( " +1) = g (n) + £g (n) until convergence is achieved. Note that in this 
relaxation scheme the potential is never used for intermediate solutions, only the potential 
increments being required. The final acceleration field is a potential gradient because the 
guess and all the intermediate solutions keep the irrotational form. The final potential field is 
computed (when needed for outputs) from the final acceleration field by numerical inversion 
of = -g. 

6. Convergence is achieved when ||5g/g||oo < s, where 5g = \\Sg\\, s is a tolerance parameter, 
and 1 1 ... i |oo is the maximum norm over all grid points. 

4. Particle-mesh scheme and time integration 

For a given set of N point particles with mass m = M/N, the Cartesian positions of the particles 
(x p ,y p ,z p ), p = 1, ...,A^, are first converted into spherical coordinates: 

r p = ^jxj+yj+zj, t P = ten-\r p /L) l/a , & p = co^\z p lr p ), <p p = tnn\y p l x p ), (25) 

where the coordinate £ is defined by equation (01. The mass of the particles is deposited on 
the radial grid, using linear or quadratic shape functions S(u - u p ) of compact support for each 
coordinate u. The resulting mass density at the grid point of indices (z, j, k) is 

M i>jJc = [pr 2 r' sin 0Af A0Ay>] . . k -rn^S (£ - £ p )S - & P )S {<p k - <p p ), (26) 

p 

where the sum extends only at the particle positions where the shape function is non zero. Since 
the mass assignment scheme is conservative, 

J] M uk = mN = M. (27) 

The inverse operation to assign the grid defined acceleration components to each particle is per- 
formed along similar lines, where now the same linear or quadratic shape functions act as inter- 
polating functions. To optimize momentum conservation, we first compute the components of 
the derivatives of the potential g = (g r , rg§, r sin §g^,) and then we interpolate them at the particle 
position x p : 

g(x p ) = Yj SUkS (ft - &)S (&j - & P )S {n - <p p ). (28) 
Finally, the interpolated spherical components 

[g r ] P = [gr] P , [g{A P = [g&/r] p , [g^p = [g^/rsin-d-],, (29) 

are combined to get the corresponding Cartesian components of the particle acceleration 

gx = [gR cos <p-gy sin <p] p , g y = [g R sin <p + g? cos tp] p , g z = [g r cos § - g# sin §] p , (30) 

where g R = g r sin § + g$ cos 

In N-MODY time integration is performed with a either second-order or fourth-order stan- 
dard leapfrog scheme. The second-order algorithm is made of the following four steps: 
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Fig. 1. CPU time (summed over all PEs) per particle per ti mestep as a function of the number of PEs for 
N-MODY simulations of MOND (red) and Newtonian (blue) Hernquist 1 199fj) spheres. 



1. half-step position move: x p (t + At/2) = x p (t) + At\ p (t)/2; 

2. evaluation of the particles acceleration at the current time: g p {t + At/2); 

3. one-step velocity move \ p (t + At) = \ p (t) + Atg p (t + At/2); 

4. a second half-step position move: x p (t + At) = x p (t + At/2) + At \ p (t + At)/2. 

The implemented fourth order leapfrog scheme is obtained by cycling three times this basic 
scheme with timestep Af replaced by (c\, Ci, c\)At respectively, where c\ = 1/(2 - 2 1 ^ 3 ) and C2 
I - 2c\. The timestep, which is the same for all particles, is adaptive in time, being determined 
by the leapfrog stability threshold Af = tj/ -y/max |V ■ g|, where the maximum is evaluated over 
all grid points and 77 is a dimensionless parameter (typically r\ = 0.3). 



5. Performances and applications 

In ICiotti et all ((2006) we tested the potential solver of N-MODY by comparing the numerical 
results with analytic MOND po tentials of sp e cial ax isymmetric and triaxial density distributions. 
The N-body code was tested in Nipo ti et all (l2007al) . where we also discussed in detail the con- 
servation of linear momentum, angular momentum and energy in N-MODY. Here we give a brief 
account of the performances of the code in typical applications. 

By choosing a tolerance parameter s = 10~ 3 for convergence in maximum norm (correspond- 
ing to =2 1CT 4 in r.m.s. norm ||...||rmsX with a relaxation parameter aT 1 = 0.3 - 0.5, the N-MODY 
solver provides the required solution in 5 - 10 iterations. For typical N-body systems, these error 
bounds correspond to an approximation in maximum norm of equation ([8]), of ||M(0)IL - 0.1 
and ||M(0)|| rms ^ 0.01. The computational time needed in the MOND potential solver scales 
roughly with the total number of grid points N g = (A^ r + l)N#N<p. In fact, in the serial version 
of the code (only one PE used), the MOND solver needs ~ 10% of the time spent to run parti- 
cles, in a typical configuration where ~ !0N g . For large grids, and lower number of particles 
per cell, some advantages are obtained by parallelization of the MOND solver. To that purpose, 
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when running with Npe processors, N-MODY adopts a domain decomposition during the iter- 
ation cycle, by assigning at each PE only a sector of the i? coordinate containing A^/TVpe grid 
points (N# > 4-Npe)- This simple strategy results to be effective, even if several operations still 
require the full grid, resulting in ~ 70% of parallelization rate. In Fig. [5]we plot the CPU time 
(summed over all PE s) per timeste p per particle as a function of the number of PEs for N-MODY 
simulations of a Hernquist ( 1990 ) sphere in Newton ian gravity (blue) and in MOND (red; with 
acceleration parameter k = 1, see Nipoti e t al. 2007c). Both simulations, using N = 10 7 particles 
and N g - 128 x 64 x 128 grid points, were run on an IBM Linux Cluster (Pentium IV/3 GHz 
PCs). In the diagram we distinguish the total CPU time (solid lines) and the CPU time spent by 
the potential solver (dotted lines). In all cases the total CPU time per particle is < 10~ 5 s, but it is 
apparent that the parallelization is not as efficient in MOND as in the case of Newtonian gravity. 
This is due to the contribution of the potential solver, which is only partially parallelized: while 
in Newtonian simulations the computational cost of the solver is always negligible, in MOND 
simulations it is non-negligible when Npe > 10, because the iterative procedure requires a factor 

of 5-10 more time than the inversion of the Poisson equation. 

We have a lready applied N- MODY to study d issipationless collap se dNipoti et al.ll2007ab . 
phase mixing dCiotti et al.l l 2007l) . galaxy merging dNipoti et al.ll2007a) and dynamical friction 
dNipoti et al.ll2008l) in MOND. We found that these dynamical processes are profoundly different 
in MOND and in Newtonian gravity. In particular violent relaxation, phase mixing and merging 
take significantly longer in MOND than in Newtonian gravity, while dynamical friction is more 
effective in a MOND system than in an equivalent Newtonian system with dark matter. 

N-MODY is publicly available upon request to the authors. It is written in FORTRAN 90 and 
can be compiled and run as either a parallel or a serial code. 

Acknowledgements. The code has been partly developed and tested at CINECA, Bologna, with CPU time 
assigned under the INAF-CINECA agreements 2006/2007 and 2007/2008. 
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